/***************************************************************************
 *   Copyright (C) 2006 by Harish Kukreja                                  *
 *   harish.kukreja@gmail.com                                              *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

//
// C++ Implementation: physiksvector
//
// Description: Definition of PhysiksBase::Vector.
//
//
// Author: Harish Kukreja <harish.kukreja@gmail.com>, (C) 2006
//
// Copyright: See COPYING file that comes with this distribution
//
//

#include <physiksbase/physiksvector.h>


PhysiksBase::Vector::Vector(unsigned int size, double init) : dimension(size) {
    cont = gsl_vector_alloc(size);
    gsl_vector_set_all(cont, init);
}

PhysiksBase::Vector::Vector(const gsl_vector *init) : dimension(init->size) {
    cont = gsl_vector_alloc(init->size);
    gsl_vector_memcpy(cont, init);
}

PhysiksBase::Vector::Vector(const Vector &init) : dimension(init.dimension) {
    cont = gsl_vector_alloc(init.dimension);
    gsl_vector_memcpy(cont, init.cont);
}

PhysiksBase::Vector::~Vector() {
    gsl_vector_free(cont);
}

void PhysiksBase::Vector::scale(double fact) {
    gsl_vector_scale(cont, fact);
}

void PhysiksBase::Vector::set_zero() {
    gsl_vector_set_zero(cont);
}

void PhysiksBase::Vector::set(unsigned int index, double value) {
    gsl_vector_set(cont, index, value);
}

void PhysiksBase::Vector::set(const Vector &copy) {
    gsl_vector_memcpy(cont, copy.cont);
}

double PhysiksBase::Vector::norm() const {
    return gsl_blas_dnrm2(cont);
}

PhysiksBase::Vector PhysiksBase::Vector::normalize() const {
    Vector res(*this);
    gsl_vector_scale(res.cont, 1.0 / gsl_blas_dnrm2(res.cont));
    return res;
}

PhysiksBase::Vector PhysiksBase::Vector::duplicate() const {
    return Vector(*this);
}

unsigned int PhysiksBase::Vector::dim() const {
    return dimension;
}

double PhysiksBase::Vector::get(unsigned int index) const {
    if (index < cont->size && index >= 0) return gsl_vector_get(cont, index);
    else return 0;
}

PhysiksBase::Vector PhysiksBase::Vector::cross(const PhysiksBase::Vector &rhs) const {
    PhysiksBase::Vector res(3, 0);
    res.set(0, this->get(1) * rhs[2] - this->get(2) * rhs[1]);
    res.set(1, this->get(2) * rhs[0] - this->get(0) * rhs[2]);
    res.set(2, this->get(0) * rhs[1] - this->get(1) * rhs[0]);
    return res;
}

PhysiksBase::Vector PhysiksBase::operator+(const PhysiksBase::Vector &lhs, const PhysiksBase::Vector &rhs) {
    PhysiksBase::Vector res(lhs);
    gsl_vector_add(res.cont, rhs.cont);
    return res;
}

void PhysiksBase::operator+=(PhysiksBase::Vector &lhs, const PhysiksBase::Vector &rhs) {
    gsl_vector_add(lhs.cont, rhs.cont);
}

void PhysiksBase::operator-=(Vector &lhs, const Vector &rhs) {
    gsl_vector_sub(lhs.cont, rhs.cont);
}

PhysiksBase::Vector PhysiksBase::operator-(const Vector &lhs, const Vector &rhs) {
    Vector res(lhs);
    gsl_vector_sub(res.cont, rhs.cont);
    return res;
}

PhysiksBase::Vector PhysiksBase::operator-(const Vector &rhs) {
    Vector res(rhs);
    gsl_vector_scale(res.cont, -1.0);
    return res;
}

double PhysiksBase::operator*(const Vector &lhs, const Vector &rhs) {
    double res;
    gsl_blas_ddot(lhs.cont, rhs.cont, &res);
    return res;
}

PhysiksBase::Vector PhysiksBase::operator*(const Vector &lhs, double fac) {
    Vector res(lhs);
    gsl_vector_scale(res.cont, fac);
    return res;
}

PhysiksBase::Vector PhysiksBase::operator*(double fac, const PhysiksBase::Vector &rhs) {
    PhysiksBase::Vector res(rhs);
    gsl_vector_scale(res.cont, fac);
    return res;
}

double PhysiksBase::Vector::operator[](unsigned int k) const {
    return gsl_vector_get(cont, k);
}

gsl_vector *PhysiksBase::Vector::core() const {
    return cont;
}
